library(corpcor)
library(GeneNet)
library(sna)
source("../er_v2.R")
source("../base_v2.R")
source("../lnet_00_functions.R")

## load
root<-"D:\\mydocs\\2012\\SNU-large-scale\\largenet\\DREAM4_Challenge2\\"
##gs.path<-"D:\\mydocs\\2012\\SNU-large-scale\\largenet\\DREAM4_Challenge2_Script\\INPUT\\gold_standards\\100"
gs.path<-"D:\\mydocs\\2012\\SNU-large-scale\\largenet\\DREAM4_Challenge2_Script\\INPUT\\gold_standards\\10"
##data.dirs<-paste(root, "DREAM4_InSilico_Size100", sep="")
data.dirs<-paste(root, "DREAM4_InSilico_Size10", sep="")
dirs<-dir(data.dirs)
gs.files<-dir(gs.path)
## for each case
##for(i in 1:length(dirs)){}
i<-1 ## for network 1
  gs.file<-paste(gs.path, gs.files[i], sep="\\")
  gs.data<-read.table(gs.file)
  wild.file<-paste(data.dirs, dirs[i], "insilico_size10_1_wildtype.tsv", sep="\\")
  wild.data<-read.table(wild.file, TRUE)
  ts.file<-paste(data.dirs, dirs[i], "insilico_size10_1_timeseries.tsv", sep="\\")
  ts.data<-read.table(ts.file, TRUE)
  ko.file<-paste(data.dirs, dirs[i], "insilico_size10_1_knockouts.tsv", sep="\\")
  ko.data<-read.table(ko.file, TRUE)
  

  ## dataset for analysis
  ## input: xx n by p where n: number of samples and p: number of variables (genes)
  xx<-ts.data[,-1]
  nvar<-ncol(xx)
  nsample<-nrow(xx)
  gene.names<-colnames(xx)

  ## True network matrix
  true.net<-matrix(0, ncol=nvar, nrow=nvar)
  colnames(true.net)<-gene.names
  rownames(true.net)<-gene.names
  for(i in 1:length(gene.names)){
    true.net[as.character(gs.data[i,1]), as.character(gs.data[i,2])]<-1
  }


## convert factor to character
i <- sapply(gs.data, is.factor)
gs.data[i] <- lapply(gs.data[i], as.character)
true.edges<-gs.data

# 
# ## output for DREAM format
# true.edges.from<-c()
# true.edges.to<-c()
# true.edges.val<-c()
# cnt<-1
# ## row
# for(i in 1:nvar){
#   ## column
#   for(j in 1:nvar){
#     if(true.net[i,j]!=0){
#       true.edges.from<-c(true.edges.from, paste("G",i, sep=""))
#       true.edges.to<-c(true.edges.to, paste("G", j, sep=""))
#       true.edges.val<-c(true.edges.val, true.net[i,j])
#       cnt<-cnt+1
#     }
#   }
# }
# true.edges<-data.frame(true.edges.from, true.edges.to, true.edges.val)



##====================== ggm =============================
## >>>>> ggm  no edges
ggm.alpha<-0.1
t.ggm1<-system.time(estimated.pcor <- ggm.estimate.pcor(xx, verbose=F))
pcor<-estimated.pcor[upper.tri(estimated.pcor)]

## edges
##ggm.edges<-matrix(0, nrow=100, ncol=3)
ggm.edges.from<-c()
ggm.edges.to<-c()
ggm.edges.pcor<-c()
cnt<-1
## row
for(i in 1:(nvar-1)){
  ## column
  for(j in (i+1):nvar){
    ggm.edges.from<-c(ggm.edges.from, paste("G",i, sep=""))
    ggm.edges.to<-c(ggm.edges.to, paste("G", j, sep=""))
    ggm.edges.pcor<-c(ggm.edges.pcor, pcor[cnt])
    cnt<-cnt+1
  }
}
## pcor test
ggm.statistic <- ggm.edges.pcor*sqrt((nsample-2-(nvar-2))/(1-ggm.edges.pcor^2))
ggm.pvalue <- 2*pnorm(-abs(ggm.statistic))
ggm.fdr<-p.adjust(ggm.pvalue, method='fdr')
ggm.edges<-data.frame(from=ggm.edges.from, to=ggm.edges.to, pval=round(ggm.pvalue,4))
ggm.edges<-ggm.edges[ggm.pvalue<ggm.alpha,]
o<-order(ggm.edges[,3])
ggm.edges<-ggm.edges[o,]
##ggm.edges[,3]<-1
write.table(ggm.edges,file="ggm_result.txt", row.names=FALSE, col.names=FALSE, quote=FALSE, sep="\t")
##time
##ggm.time<-t.ggm1[1]+t.ggm2[1]+t.ggm3[1]


## convert factor to character
i <- sapply(ggm.edges, is.factor)
ggm.edges[i] <- lapply(ggm.edges[i], as.character)
test.edges<-ggm.edges
##test.edges<-true.edges[1:13,]

## get
ggm.evals<-net.evaluation(true.edges, test.edges)




##------------enet-----------
selbs<-matrix(0, nvar,nvar)
colnames(selbs)<-gene.names
rownames(selbs)<-gene.names
for(i in 1:p){
  y<-x[,i]
  if(var(y)!=0){
    dat<-x[,-i]
    fit<-glmnet(dat, y, lambda=lam, alpha=alpha)
    selbs[-i,i]<-as.numeric(fit$beta)
    ##cat(sum(selbs[,i]!=0),"---", i, "\n");flush.console()
  }
}


 
